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Summary 


This report describes the assignment posed by NASA/Goddard Space Flight Center. The 
assignment had three parts: 

• Comparing particle trajectories computed, first in an analytically specified set of global 
electric and magnetic fields; and then, in the same fields as sampled onto a finite element 
mesh in three dimensions, with interpolation of intermediate values. 

• Studying the effect of mesh spacing, resulting in an evaluation of adequate spacing 
resolution. 

• Studying time-dependent fields in the context of substorm dipolarizations of the 
magnetospheric tail. 

One of the aims of the assignment was to gain confidence when it comes to executing MHD- 
simulations, by investigating the interpolating simulations. Interpolations are much faster to do 
than analytical simulations. The study was designed to validate and document work that NASA 
has been doing to simulate test particle motions in complex three-dimensional 
magnetohydrodynamic (MHD) simulations of the Earth’s magnetosphere. 

During the whole process, the test particle code of Delcourt was used. The three-dimensional test 
particle trajectory code of Delcourt uses the full equation of motion, and thus, including drift 
velocities of second order, in order to examine in detail the motion of near-Earth (~10- 15 Re) 
plasma sheet particles during substorms. The emphasis is placed on orbit features in the course of 
expansion phase, thus ignoring the subsequent evolution of the ion distributions produced. A total 
of 18 different scenarios were used during the different phases of the work. The scenarios 
differed in all factors, except the Tsyganenko level and the dipole tilt angle, which were kept 
constant except for the dynamic simulations where the Tsyganenko levels were also varied. 

The results from the particle motion in an interpolated magnetic field and particle motion in 
interpolated magnetic and electric fields showed a very good resemblance and thus, we can be 
confident when it comes to executing MHD- simulations (with respect to interpolating 
simulations). Thus, a lot of time can be saved by using interpolation — instead of time consuming 
analytical simulations — in these cases. The sensitivity analysis — for the particle motion in an 
interpolated magnetic field and particle motions in interpolated magnetic and electric fields — in 
general show large improvements in accuracy between the analytical results and the precalculated 
results with an increased number of grid points. In some cases, more or less identical results 
between the analytical and the precalculated simulations were found. The decreased accuracy 
concerning energy could be explained by the random nature of the particle trajectories. Weighing 
the simulation time (the simulations take more time when increasing the number of grids) and the 
resulting accuracy, an optimal number of grids would probably be eight times the original 
number. 

The few things that differed when comparing the polynomial results with the linear results — 
dynamic magnetic and electric fields — were generally easily corrected by increasing the time 
interval between the change of the Tsyganenko levels. By doing so, the time interval will make 
the curve much smoother and the transition less abrupt. The comparison concerning a dynamic 
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magnetic field showed a very good resemblance between the two cases, but when trying to add a 
dynamic electric field, the mathematics posed a problem. The comparison did not show enough 
resemblance between the two cases. The problem was finally solved, and when comparing the 
results of a dynamic magnetic and electric field, a very good resemblance with the analytical 
results was found, thus, satisfactory results were achieved. 
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1. BACKGROUND 


1.1 Assignment 

The assignment was performed at NASA/Goddard Space Flight in Greenbelt, Maryland. The 
assignment had three parts: 

• Comparing particle trajectories computed, first in an analytically specified set of global 
electric and magnetic fields; and then in the same fields as sampled onto a finite element 
mesh in three dimensions, with interpolation of intermediate values. 

• Studying the effect of mesh spacing, resulting in an evaluation of adequate spacing 
resolution. 

• Studying time-dependent fields in the context of substorm dipolarizations of the 
magnetospheric tail. 

One of the aims of the assignment was to gain confidence when it comes to executing MHD- 
simulations, by investigating the interpolating simulations. Interpolations are much faster to do 
than analytical simulations. The study was designed to validate and document work that NASA 
has been doing to simulate test particle motions in complex three-dimensional 
magnetohydrodynamic simulations of the Earth’s magnetosphere. 


1.2 Theory 


1.2.1 Ionosphere 

The ionosphere is a neutral atmosphere that contains a source of ionization for the gases in its 
atmosphere. The sources of ionization include photons and energetic-particle “precipitation.” The 
photons come primarily from the Sun as EUV/UV radiation. The process involving photons is 
called “photoionization” and the process involving energetic particles is called “impact 
ionization.” Ionizing particles may come from the galaxy (cosmic rays), the Sun, the 
magnetosphere, or the ionosphere itself if a process for local ion or electron acceleration is 
operative. Precipitating energetic electrons produce additional ionizing photons within the 
atmosphere by a process known as Bremsstrahlung, or braking radiation. The ionization due to 
particles causing aurora, takes place mainly in the altitudes from 90-300 km. Cosmic rays with 
energies in the mega electron volt and giga electron volt intervals, cause ionization at altitudes 
from 90-70 km. [2] 

The ionosphere consists of three major layers or regions: the D region (below 90 km), the E 
region (from 90-130 km), and the F region (above 130 km). The F region is usually further 
divided into FI and F2 layers, because a second ledge sometimes appears in its profile below the 
main (F2) peak. One can think of the layers as being independently produced by the absorption of 
solar radiation by specific components of the neutral atmosphere, which respond differently to 
different parts of the incoming solar photon spectrum. The E layer is usually clearly noticeable as 
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a change in slope in daytime electron density profiles near 110 km. The ions in this layer are 
mainly O 2 and NO + , which have been produced by UV radiation in the 100-150 nm range and 
solar x-rays in the 1-10 nm range. The peak density of this layer is close to the peak of the 
production rate for these ions. Vertical transport of ions is considered to play a minor role in the 
formation of this layer. The F 1 layer is composed primarily of 0 + . The maximum electron density 
in this layer occurs at approximately 170 km, which is close to the level of the maximum ion 
production by photons in the spectral range from about 17-91 nm. The F2 layer contains the main 
peak in the ionosphere density. This peak, due to collisions at lower altitudes, keeps the density 
of ionized particles down, while the density at higher altitudes is kept down by the lack of 
particles to be ionized. [2] This peak density is also in a region dominated by O . The lowest 
ionosphere, the D region, is of practical interest because of its role in commercial radio 
communication. The high ion-neutral collision frequencies make radio-wave absorption there 
important, and so the electron density is of primary concern. Only the most energetic ionization 
sources can penetrate to D-region altitudes. Between about 80-90 km, 0.1-1 nm x-rays from the 
Sun are the primary sources; the very intense Lyman-a (121.6 nm) radiation from the Sun has its 
peak ion production rate at about 70-80 km, and ionization of cosmic-ray particles dominates 
below. The predominant ions, NO 4 and O 2 , can recombine with electrons, but at these low 
altitudes the electrons can also attach themselves to neutrals to form negative ions. In general, the 
variation of ionospheric density during the day depends on altitude through both the local ion 
composition and the source variation during the day. [1] 
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Figure 1. Ionospheric and atmospheric structure. [7] 


1.2.2 Magnetosphere 

The magnetosphere consists of the magnetopause, magnetosheath, and the magnetotail. The 
magnetopause and the magnetotail are the two regions within the Earth’s magnetosphere where 
relatively thin sheets of electric current separate regions of different magnetic fields. The 
magnetopause is usually placed at a distance of 10 Earth radii on the dayside. The magnetosheath 
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is found within the bow shock. In the magnetosheath, the solar wind plasma flows with a 
subsonic speed. The magnetotail is found on the nightside. [2] 
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Figure 2. The Earth’s magnetosphere. [7] 

A current sheet can be defined as a thin surface across which the magnetic-field strength and/or 
direction can change substantially. Current sheets come naturally out of the frozen- in- flux 
concept. Collisionless plasmas do not mix easily; instead, they tend to form cells of relatively 
uniform plasma having a magnetic field flowing through it. These cells are separated by thin 
current sheets through which little or no magnetic flux crosses. The interaction between different 
plasma regimes occurs at these thin boundaries. Magnetic reconnection is one of these processes; 
see below for magnetic reconnection. 

The geomagnetic tail is the region of the Earth’s magnetosphere that stretches away from the Sun 
behind the Earth. It acts as a reservoir of plasma and energy. The energy and the plasma are 
released into the inner magnetosphere aperiodically during magnetospheric substorms. A current 
sheet is situated in the center of the tail, embedded within a region of hot plasma — the plasma 
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sheet — which separates two regions called the tail lobes. The two tail lobes connect magnetically 
to the two polar regions of the Earth and are identified as the north and south lobes. 

The highly stretched magnetic-field geometry of the Earth’s magnetotail implies a westward tail 
current across the center of the tail, near the equatorial plane. The inner magnetosphere contains a 
ring current, which flows westward in approximate circles centered on the Earth. Additional 
partial-ring currents also flow partway around the Earth in the middle magnetosphere. Birkeland 
currents connect the ends of the partial rings to the ionosphere, where conduction currents can 
complete the circuit. Additional Birkeland currents connect to the tail and magnetopause 
boundary layer. [1] 

The magnetopause is the upper boundary of the magnetosphere. It separates the geomagnetic 
field and plasma of primarily terrestrial origin from solar-wind plasma. In the simplest 
approximation, the magnetopause can be considered as a boundary separating a vacuum magnetic 
field from plasma. Three distinct types of magnetopause boundary layers are recognized: 

• The Low-Latitude Boundary Layer (LLBL), a region that contains a mix of 
magnetosheath and magnetosphere plasma, and within which plasma flows can be found 
in almost any direction, but are generally intermediate between the magnetosheath flow 
and magnetospheric flows. The entry layer is found in the region where the magnetic null 
or cusp occurs in a closed magnetopause model and extends about 3h either side of noon. 
The plasma is characteristic of the magnetosheath, but the flows are low-speed and 
disordered. 

• The High-Latitude Boundary Layer (HLBL), or plasma mantle, is found at higher 
latitudes tailward of the cusp and entry layer. Within HLBL flows are always tailward, but 
flow speed, density, and temperature all decrease away from the magnetopause. In 
general, the HLBL is spatially uniform, with a gradual transition from magnetosheath to 
lobe properties. It often has no distinct inner edge, in sharp contrast to the LLBL, which 
usually has a very distinct inner edge. The entry layer and plasma mantle are generally 
agreed to be on open magnetic-field lines. They are populated by a mixture of 
magnetosheath plasma, that entered the magnetosphere along open field lines in the cusp 
and ionospheric plasma, which flowed up from the cusp and polar cap in an upward flow 
known as the “polar wind.” Reconnection is assumed to occur at the nose of the 
magnetopause. Magnetosheath particles flow along the newly opened field lines. The 
particles mirror as they encounter the stronger magnetic field nearer the Earth. After 
mirroring, they move back up the field line, joined by lower-energy ionospheric particles. 
The field line has meanwhile convected tailward to become a lobe field line, still near the 
magnetopause, or a mantle field line. Lower-energy particles move slower and thus, take 
longer to reach a given distance downtail. In this longer time, the field line will have 
convected farther from the magnetopause. 

• The Plasma-Sheet Boundary Layer (PSBL) forms the boundaries between the plasma 
sheet and the two tail lobes. Within this boundary layer, plasma-number density and 
temperature are intermediate between lobe and central plasma sheet values, but its most 
characteristic feature is sustained field-aligned ion and electron flows, directed both 
Earthward and tailward. The PSBL is formed of plasma that has just been accelerated in a 
reconnection region. 


4 



The most important currents that give rise to the magnetospheric magnetic field are those that 
flow inside the Earth. The Earth’s field is effectively a perfect dipole at distances more than about 
two Earth radii from the center of the Earth. Plasma flows over the polar-regions from noon 
toward midnight and then back toward the dayside at somewhat lower latitudes, corresponding 
roughly to the auroral zone, on both the dawn and dusk sides (=“convection”). 
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1.2.3 Magnetic Reconnection 


Magnetic reconnection is a process that can take place only if the magnetic-field lines being 
convected with the fluid breaks down. Magnetic-field lines enter the diffusion region from the top 
and bottom, and instead of being annihilated, they leave from both sides. In the process, they are 
cut and reconnected to different partners. Plasma that was originally on different flux tubes, 
coming from different regions, now finds itself on a single flux tube in total violation of the 
frozen-in-flux theorem. Previously, the current sheet separated two magnetic-field regions; now 
magnetic flux crosses the current sheet. In other words, we now have an open boundary with a 
finite normal magnetic-field component, as opposed to the previous closed boundary. Plasma can 
cross the boundary by following flux tubes. 


• Magnetic field. H 



plasma flow 


Figure 3. Magnetic reconnection at the diffusion region. [8] 


1.2.4 Ionosphere-Magnetosphere Coupling 

The role of the ionosphere as a plasma source for the magnetosphere has been increasingly 
considered since the discovery of energetic 0 + in the auroral zone. Observations throughout the 
magnetosphere have revealed the presence of ionospheric ions in all regions: auroral zone, 
plasmasphere, plasma sheet and lobe, magnetosheath, and boundary layers. [6] The key 
electrodynamic element in the physical coupling between the ionosphere and the magnetosphere 
is the current that flows along field lines between the two regions. 

For typical steady conditions, the inner edge of the plasma sheet can effectively shield the inner 
magnetosphere, and the low-latitude and mid-latitude ionosphere to which it is connected, from 
the convection electric field. Shielding by the inner edge of the plasma sheet is frequently 
ineffective. Time variations in the applied convection potential or other magnetospheric 
parameters can cause temporary penetration of the shielding. It is also possible that the plasma 
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sheet is sometimes incapable of effective shielding for longer periods. Furthermore, neutral winds 
generated by magnetospheric activity can cause significant electric fields at low and middle 
latitudes in the ionosphere. Significant potential drops are observed to occur along auroral-zone 
magnetic-field lines. These potential drops usually occur in regions of strong upward field- 
aligned current. 

Two major processes cause magnetospheric particles to be lost in the Earth’s atmosphere, 
namely, pitch-angle, scattering, and charge exchange. A particle that conserves its first adiabatic 
invariant is likely to be lost to the atmosphere if its mirror magnetic field exceeds the magnetic 
field at a few hundred kilometers altitude. Auroral-particle precipitation proceeds primarily 
because plasma waves change particle pitch angles in violation of the first adiabatic invariant. 
Pitch-angle scattering occurs much more slowly in the radiation belts than in the plasma sheet, 
and radiation-belt particles can orbit the Earth for days, weeks, or months, depending on their 
species, energy, and location. [1] 

Cold plasmas from the ionosphere contribute to the hot plasmas of the magnetosphere. On openly 
convecting field lines, polar wind occurs continuously as convection opens the field lines and 
empties their accumulations into the polar lobes and downstream solar wind, so that they never 
reach equilibrium pressures and are therefore called “refilling flows.” The outer plasmasphere has 
been shown to flow Sunward during magnetospheric disturbances and these cold plasmas have 
been discovered to be present in the subsolar low-latitude magnetopause region under a wide 
variety of conditions. When strong convection drains away part of the plasmasphere, the supply 
of plasma is enhanced in a transient way by the rapid release of accumulated plasma. Under 
steady conditions, however, the plasmasphere remains trapped and the magnetosphere is supplied 
only from the higher-latitude polar wind regions. The contribution of ionospheric light ions to 
magnetospheric hot plasmas is not as well established and is complicated by the difficulty of 
discriminating protons of solar or geogenic origin. 

The light ion auroral outflows have fluxes similar to polar wind, while the heavy ion outflows 
have fluxes that range from much less than, to much greater than, polar wind outflows depending 
on the free energy available. [3] Heavy ions are important during magnetospheric storms, when 
the dissipation of energy in the ionosphere proper is substantial. Part of this energy goes into 
energization of heavy ions sufficient to overcome their gravitational binding to the Earth. Other 
research has shown that the supply of ionospheric plasma to the plasma sheet, especially heavy 
0 + plasma, is not only important, but is strongly modulated by convection as driven by the 
Interplanetary Magnetic Field. [3] 


1.2.5 Adiabatic Invariants 


If the changes in the magnetic field encountered by a particle within a single gyration orbit, will 
be small compared with the initial field, then the particle’s magnetic moment, p, is given by 


h = 


mv\ 

2 Tfi ’ 


( 1 ) 
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where m is the mass of the particle, v ± is the velocity perpendicular to the magnetic field, and the 
quantity B is the magnetic field. The magnetic moment will remain constant, and the quantity u is 
also called the first adiabatic invariant. Adiabatic in this case refers to the requirement that u may 
not remain invariant or unchanged unless the parameters of the system — such as its field strength 
and direction — change slowly. The second adiabatic invariant is also called the bounce invariant, 
which is given by 


J = 


ni\ 


j ) p dl = 2 jmv^dl , 


( 2 ) 


where p« is the component of the momentum along the magnetic field, v,, is the component of 

the velocity along the magnetic field, / is the distance along the field line, and mi and nij are the 
locations of the particle’s mirror points. 

What this means is that a particle trapped on magnetic field lines and conserving its adiabatic 
invariants will be confined to a surface specified by the distance / to the equatorial crossing of the 
field line over half a bounce from mirror point mi to nij. 

As a particle bounces and drifts in a magnetic field that varies slowly along its orbit, it traces out 
a so-called “drift shell.” The flux that is enclosed by this shell is the same for any surface that 
intersects along a closed curve. If the field changes slowly (on time scales long compared with 
the time required for the particle to drift fully around the field), the flux enclosed by the drift shell 
will remain unchanged. The third adiabatic invariant is defined as the magnetic flux enclosed by 
the drift shell of a particle. [1] 


1 .2.6 Magnetohydrodynamics (MHD) Theory 

The relevant fluid equations for plasmas are the equations of hydrodynamics, but because electric 
and magnetic fields and currents are always important in plasmas, one must include their effects 
as well. This means that the equations of magnetohydrodynamics (MHD) must be introduced. 
The equations incorporate familiar mechanical laws and also account for electromagnetic 
properties. 


The continuity equation, which states what happens to the number density of the fluid as it moves 
from one place to another, says that the number of particles can change only if there are particle 
sources ( S s ) or losses (L s ) which are given by 


dry 

dt 


+ V • n u = S - L, 


( 3 ) 


where n s is the number density, t is the time, u s is the flow velocity, S s is the particle sources, and 
L s is the particle losses. 

The momentum equation is given by 
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( rj(P ; Uv) + V • (p v u v u A )) = -Wp s + p E + j s x B + p s F / m s , (4) 
or 


where p (/s is the charge density, j s is the current density, p s is the mass density, p s is the pressure 
force, E is the electric field, B is the magnetical field, F„ is the non- electromagnetical forces 
(such as gravitational forces), and m s is the particle mass. 

In order to consider the equations that describe how current and charge density affect the 
magnetic and electric fields, Maxwell equations are used. Included in the Maxwell equations is 
Poisson’s equation — which deals with the electric field — which is given by 

V-E = P,/ £o , (5) 


where p q is the net charge density and So is the permittivity of free space. 


Faraday’s law — which deals with the electric field — is given by 


SB 

dt 


= -V x E . 


( 6 ) 


Ampere’s law — which relates the magnetic field to the net current — is given by 


dE 

V x B = ,u 0 ( j + £ 0 — — ) , 
at 


( 7 ) 


where p 0 is the permeability of free space and j is the net current. 
The requirement that B is divergenceless, given by 


V -B = 0 . (8) 

Lorentz-force law is given by 

F = g(E + uxB) . (9) 

When E + uxB = 0 is valid, one can demonstrate that in MHD fluid flow, the magnetic flux can 
be frozen into the fluid. The frozen- in-flux conditions tells us that if we follow the fluid initially 
on a certain surface as it moves through the system, the flux through the surface will remain 
constant even as the surface changes its location and its shape. [1] 
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Figure 4. MHD simulation of Earth’s magnetosphere. [9] 
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2. Method and Results 


2.1 Computer Tools and Codes Used 

During the work, the LINUX operating system was used. The simulations and calculations were 
done using FORTRAN-codes. The graphical presentations were made by using Interactive Data 
Languages (IDL) code, which is generally used by NASA. 


2.1.1 FORTRAN 

The name FORTRAN is derived from FORmula TRANSlation, indicating that the language was 
intended from the start for translating scientific equations into computer code. The first version of 
the FORTRAN language was developed during the years 1954-1957 by IBM. When FORTRAN 
was introduced, it was a success from the start as it made programming so much easier than 
machine language. FORTRAN has gone through changes throughout the years, being updated in 
1962, 1977, 1990, 1996, and at present (FORTRAN 2000). [4] FORTRAN is still indeed a living 
language and still very much in use in the scientific world. One of the reasons why FORTRAN is 
still in use is that so much scientific material is written in FORTRAN and to translate this into 
another language would simply cost too much. 


2.1.2 IDL 

IDL stands for Interactive Data Language, which is an array-oriented data analysis and 
visualization environment. The first version of IDL was released in 1981. Since then, IDL has 
emerged from its roots in the astronomical and space sciences to become a widely used tool in 
research, educational, and commercial settings in areas as diverse as the Earth sciences, medicine, 
physics, and engineering tests and analysis. IDL is used, for example, when analyzing and 
visualizing images from the Hubble Space Telescope. One of the strengths of IDL is its support 
for a variety of hardware and operating systems combinations. [5] 


2.1.3 Test Particle Code ofDelcourt 

During the entire process, the test particle code of Delcourt was used. The three-dimensional test 
particle trajectory code of Delcourt uses the full equation of motion and thus it includes drift 
velocities of the second order, in order to examine in detail the motion of near-Earth ( — 10 — 15 Re) 
plasma sheet particles during substorms. The emphasis is placed on orbit features in the course of 
the expansion phase. 

In the test particle code ofDelcourt, the magnetic field model of Tsyganenko 89 was used. The 
magnetic field B was therefore considered to be the sum of the Earth dipole, B dip , and an external 
contribution, B ext . In order to include the electric field induced by the time-varying magnetic 
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field, a vector potential A associated with the perturbation B ext was calculated. The instant 
magnetic field and induced electric field in a position r are given by 

B(r,t) = B dip (r) + curl[A(r,t)] (10) 

and 

E ind (r,t) = -dA(r,t)/dt . (11) 

The full equation of motion was used in order to track all ion drift paths and was given by 

r = (q/m) [E + r xBj +g , (12) 

where r is the instant position of the particle, q is the charge of the particle, m is the mass of the 
particle, g is the gravitation, and E is the combined induced and polarization electric fields. 

Equation (12) was integrated using the Runge-Kutta technique of the fourth order. The above 
described set of equations, was solved in the FORTRAN-code used during all simulations (with 
conditional conditions being fulfilled). 

During the expansion phase of a substorm, single-particle trajectory simulations demonstrate 
increased Earthward convection of low-latitude plasma in the midtail region (~10 Re). They 
indicate possible dramatic ion accelerations, mainly in the perpendicular direction. They 
underline the determining role of the particle drifts, which have usually been considered 
negligible and show clearly the following two effects: 

1. A possible creation of new high-altitude mirror points results from large centrifugal 
accelerations across the magnetospheric field lines and is controlled by particle parallel velocity. 

2. A possible breaking of the first adiabatic invariant relates to intense polarization drifts and 
depends on particle mass and/or charge state via gyroperiod. 

Systematic orbit calculations show that the low-latitude populations are most affected by the 
convection of the geomagnetic tail and reveal possible dramatic accelerations of the particles. 
These trajectory results are of relevance to explain various storm time signatures at 
geosynchronous altitudes and provide a populating mechanism for the ring current region. [11] 


2.2 Particle Motions in Interpolated Magnetic Field 

The first simulations were aimed at getting the results of the analytical cases to fit the results of 
the precalculated cases, with only a magnetic-field added (only co-rotation electric field was 
included) and with no variations in time. Tsyganenko level 2 in the Tsyganenko 89 model 
(describing geomagnetic conditions) was used. Tsyganenko ’s magnetospheric magnetic field 
model of 1989, takes into account the effect of warping the tail current sheet in two dimensions 
due to the geodipole tilt, as well as spatial variations of the current sheet thickness along the Sun- 
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Earth and dawn-dusk directions. The corresponding analytic forms for the magnetic field 
components are obtained using an indirect approach in a two-stage procedure. [12] The tilt angle 
of the dipole was set to 0°, the simulated particle was a proton, the initial particle energy was set 
to 1000 eV, and the initial Magnetic Local Time (MLT) was on the nightside during all 
simulations. MLT at a given location is determined by the angle subtended at the geomagnetic 
axis between the geomagnetic midnight meridian and the meridian that passes through the 
location. The following factors were varied: initial distance to Earth, initial magnetic latitude, 
initial pitch angle, and initial gyrophase. See Tables 1 and 2, for a full description of the 
simulated scenarios that were used throughout the report. 
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In order to get the results of the analytical cases to fit the results of the precalculated cases, the 
following measures were taken: 


1 . Only the co-rotational part of the electric field was kept, the convectional part was set to zero. 

2. Only the external field was interpolated, the dipole field was then added on. 

3. A spherical grid-system was used. 

A total of 1 8 different scenarios were used during the different phases of the work. The scenarios 
differed in all included factors except the Tsyganenko level and the dipole tilt angle, which were 
kept constant except for the dynamic simulations where the Tsyganenko level was also varied. 
The Tsyganenko level describes the geomagnetic conditions using a scale from 1 to 6, where 1 
equals quiet geomagnetic conditions and 6 equals the most disturbed geomagnetic conditions. 

a. The dipole tilt angle is the angle between the geomagnetic dipole and the plane orthogonal to 
the Earth-Sun line (varies within a range from -35° to +35°). 

b. The magnetic local time will vary between 00 hours to 24 hours. 

c. The magnetic latitude will vary between -90° to +90°. 

d. The pitch angle is the angle between v and B. 

e. The pitch angle will vary between 0° and 180°. 

f. The gyrophase will vary between 0° and 360°. 

Table 1. Scenario 1 through 9 


Scenario 

1 

2 

3 

4 

5 

6 

7 

8 

9 

Tsyganenko 

level 

2 

2 

2 

2 

2 

2 

2 

2 

2 

Dipole tilt 
angle (°) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Initial 

geodistance 

(RE) 

2 

4 

6 

8 

10 

12 

14 

6 

2 

Initial 
magnetic 
local time 
(h) 

20 

20 

23 

23 

23 

23 

23 

23 

23 

Initial 
magnetic 
latitude (°) 

45 

45 

45 

45 

45 

45 

45 

10 

10 

Initial 

energy (eV) 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

Initial pitch 
angle (°) 

45 

45 

45 

45 

45 

45 

45 

45 

45 

Initial 

gyrophase 

(°) 

45 

45 

45 

45 

45 

45 

45 

45 

45 
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Table 2. Scenario 10 through 18. 


Scenario 

10 

11 

12 

13 

14 

15 

16 

17 

18 

Tsyganenko 

level 

2 

2 

2 

2 

2 

2 

2 

2 

2 

Dipole tilt 
angle (°) 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Initial 

geodistance 

(RE) 

12 

2 

6 

12 

2 

6 

12 

6 

12 

Initial 
magnetic 
local time 

(h) 

23 

23 

23 

23 

23 

23 

23 

23 

23 

Initial 
magnetic 
latitude (°) 

80 

45 

45 

45 

45 

45 

45 

10 

10 

Initial 

energy (eV) 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

2000 

3000 

Initial pitch 
angle (°) 

45 

10 

20 

80 

45 

45 

45 

45 

45 

Initial 

gyrophase 

(°) 

45 

45 

45 

45 

10 

20 

80 

45 

45 


Very good results were achieved for the simulated scenarios, where the energy in both the 
analytical and the precalculated cases were kept relatively constant and the particle trajectories 
showed very similar results, which was desirable. See Figures 5 and 6 for an example (Scenario 
4). As seen from Figures 5 and 6, the energy is kept fairly constant in both cases. The pitch angle, 
magnetic moment, and to a lesser degree, the kappa-value are very similar in both cases. The 
kappa-value is the square root of the magnetic field curvature radius to the particle gyroradius. 
[13] The appearance of the trajectories has similarities. 
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Figure 5. Results from the analytical case, interpolated magnetic field (Scenario 4). 
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Figure 6. Results from the precalculated case, interpolated magnetic field (Scenario 4). 


2.3 Particle Motions in Interpolated Magnetic and Electric Fields 

The work now continued with adding an electric field. A FORTRAN-code was developed, which 
generates an electric field in accordance to Stem-Volland. The Stem-Volland electric field 
model is K p -dependent. Later models are also time-varying with respect to the convection electric 
field. [10] During the development of the code, the region below the ionosphere was assumed to 
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be a dipole. The generated electric field was used in the simulation code, where the co-rotational 
part was also activated. Particle trajectories with initial conditions listed in Tables 1 and 2 are 
calculated in the precalculated magnetic and electric fields. The trajectories are compared with 
those with the same initial conditions, but are calculated in analytical fields. 

In order to get an idea of the procedure of the simulations, Figure 7 will give an idea of how the 
different codes were used during the simulations. 



Figure 7. The scheme of the different codes involved. 

The results of the analytical and the precalculated cases were compared. Minor adjustments 
concerning the sampling interval of the electric field were made in order to improve the 
resemblance between the cases. The results showed a very good resemblance, concerning energy, 
pitch angle, magnetic moment, kappa-value, and the appearance of the trajectory. See Figure 8 
for an example (Scenario 3). 
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Figure 8. Analytical case to the left and the precalculated case to the right (Scenario 3). 


2.4 Dynamic Magnetic and Electric Fields 

During the simulations, while varying the magnetic and electric field, a linear and a polynomial 
transition between Tsyganenko level 5 to level 2 was used. The results of the two different 
transitions were then examined in order to look for differences in the resulting parameters. As for 
the linear and polynomial transition, we use the function B(t) given by 

B(t) = B(t 0 ) + F(t)(B(t f )-B(t 0 )), (13) 


where the linear transition uses: 


F{t) = — — , (14) 

tf t 0 

while the polynomial transition — in this case — uses: 

F(t) = 10- -15- + 6 • (^-^) 5 . (15) 

tf —t 0 tj — t 0 tf — 1 0 
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After comparing the simulations with the two different transitions, the work started on developing 
a trajectory code where a dynamic magnetic field (level 5 to level 2) and interpolation 
(precalculated) was used. When the results were found to be good enough, the electric field was 
to be added on. 

When comparing the results of the linear transition and the polynomial transition, a very good 
resemblance was found (nearly identical resemblance), concerning pitch angle, kappa-value, and 
the appearance of the trajectory. The only things that differed were the magnetic moment — where 
the linear transition results in a sudden and powerful change of the magnetic moment, which 
depends on the sudden change in the electric field — and the energy to a certain amount. In one 
case (Scenario 8), the polynomial transition showed a sudden and powerful change in the kappa- 
value, which the corresponding linear transition did not show. The time interval of the transition 
was increased in this case from 60 to 600, which resulted in the powerful change of the kappa- 
value disappearing. Increasing the time interval means that the transition is being made more 
smoothly. See Appendix 1 for the linear and the polynomial transition for Scenario 8. See 
Appendix 2 for the linear and the polynomial transition for Scenario 8, when increasing the time 
interval to 600. 

When comparing the results of a dynamic magnetic field, a very good resemblance with the 
analytical results was found, except for one scenario (Scenario 10). In this scenario, the particle 
was found in the polar cusp and this is where the magnetic field has discontinuities, which will 
make interpolation hard to achieve. See Appendix 3 for graphs of Scenario 10. 

When the electric field was added, the newly generated electric field was checked against the 
earlier electric field by checking the potential in a certain point. The values of the potential in a 
point were found to be the same for both fields. See Figure 9, for a view of the electric fields in 
X-Z-plane, e x -direction. 
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Figure 9. The electric fields in X-Z-plane and e x -direction. 

Simulations were now conducted with a dynamic magnetic and electric field added on. The 
results between the analytical and the interpolated (precalculated) simulations were found to have 
good resemblances. The resemblances concerning the electric field, however, were not found to 
be adequate, thus the work continued with the electric field. The problem was how to calculate 
the expression given by 


f( E t »B)ds , (16) 

without field line tracing. Equation (16) is found in the expression for the total E-field, given by 

E = -W(P+f(E i »B)ds) + E t± , (17) 

where P is the electrostatic potential (obtained from mapping to the ionosphere), E, is the 
inductive electric field, and B is the magnetic field. The problem was solved in the following 
way: the total E-field was precalculated during the transition between the two levels, calculating 
the E-field at a 3 second interval. The particle trajectory was calculated by doing temporal and 
spatial interpolation to find the E-field at a particular time and location along the trajectory. When 
comparing the results of a dynamic magnetic and electric field, a very good resemblance with the 
analytical results was found. See Figures 10 and 1 1 for an example (Scenario 1). 
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Figure 10. Analytical case (Scenario 1). 


22 


2 r 

i ■ 

o - 


Test-particle code of D. Delcourt 


M _ 2 


-3 
-4 L 



Figure 11. Precalculated case (Scenario 1). 


3. Analysis and Discussion 


3.1 Particle Motions in Interpolated Magnetic Field 

The measures taken were found to be adequate to get the analytical cases to fit the results of the 
precalculated cases. The results compared were the trajectories, energy, pitch angle, magnetic 
moment, and the kappa-value. The most interesting input values were varied in the different 
scenarios. They were varied within a quite diverse range of numbers (see Tables 1 and 2). Still, 
the results showed very good resemblance between the analytical cases and the precalculated 
cases. 
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3.2 Particle Motions in Interpolated Magnetic and Electric Fields 

The measures taken were found to be adequate to get the analytical cases to fit the results of the 
precalculated cases. The results compared were the trajectories, energy, pitch angle, magnetic 
moment, and the kappa-value. The same accuracy in resemblance was received as for the case 
with the interpolated magnetic field. 


3.3 Sensitivity analysis 

Simulations were carried out for all scenarios — with particle motions in interpolated magnetic 
and electric fields — using grid systems, 8 times the original number of grids, 27 times, and finally 
half the number of grids. A simulation using 64 times the original number of grids was stopped as 
the computer crashed. When developing the grid systems, the given grid points from the original 
grid system were used. We obtained the new grid points using linear interpolation. See Figures 12 
and 13 for the original grid system and the grid system with 8 times the original number of grids. 


spherical Grid 



x m 


Figure 12. The original grid system. 
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Figure 13. The grid system 8 times the original number of grids. 


The following results were received for the different scenarios: 
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Scenario 1: The resemblance to the analytical case gets better when it comes to pitch angle, 
magnetic moment, and the appearance of the trajectory, but unfortunately, in the case of a grid 27 
times the original amount, the particle precipitated shortly after release and no conclusions can be 
made. 

Scenario 2: The resemblance to the analytical case gets better (closer resemblance between 
analytical result and precalculated result when increasing the amount of grids) when it comes to 
energy, kappa-value, pitch angle, and the appearance of the trajectory. 

Scenario 3: Same as for Scenario 2. See Appendix 4 for the graphs of Scenario 3. 

Scenario 4: The resemblance to the analytical case gets better for every increase in the amount of 
grids when it comes to magnetic moment. Concerning the energy and pitch angle, the result gets 
better with every increase except for the last one (27 times), when the tendency goes somewhat in 
the wrong direction. 

Scenario 5: Same as for Scenario 2, except that the energy goes somewhat in the wrong direction 
for the case when the grid amount is 27 times the original amount. 

Scenario 6: Because the particle goes out of the predefined boundary for all precalculated cases at 
an early stage, no conclusions regarding the development of the resemblance can be made. 

Scenario 7: Because the particle goes out of the predefined boundary for all cases at an early 
stage, no conclusions regarding the development of the resemblance can be made. 

Scenario 8: Same as for Scenario 5. 

Scenario 9: Same as for Scenario 2, but the particle impacts for all precalculated cases. 

Scenario 10: Same as for Scenario 7. 

Scenario 1 1 : Same as for Scenario 7, except that the particle impacts. 

Scenario 12: Same as for Scenario 2. 

Scenario 13: The resemblance decreases as the particle goes out of the predefined boundary as the 
amount of grids increases. 

Scenario 14: No greater improvement in resemblance is detected. 

Scenario 15: The results somewhat oscillate. The results of the normal grid amount show good a 
resemblance. It gets somewhat worse for the grid with 8 times the normal amount, and then gets 
better for the grid with 27 times the normal amount. 
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Scenario 16: The resemblance gets much better when increasing the amount of grids 27 times, as 
the particle goes out of the boundary for the other precalculated results. See Appendix 5 for the 
graphs of Scenario 16. 

Scenario 17: Same as for Scenario 2, except that the magnetic moment gets somewhat worse 
when increasing the amount 27 times. 

Scenario 18: Same as for Scenario 7. 


3.4 Comparing the Polynomial Results with the Linear Results (Dynamic Magnetic and 
Electric Fields) 

When comparing the polynomial results with the linear results, the results were almost identical. 
The only things that differed were the magnetic moment — where the linear result experiences an 
abrupt change, which is due to an abrupt change in the ii-field — and/or possibly, the energy. For 
Scenario 8, the polynomial result showed an abrupt change in the kappa-value, which the linear 
result did not show. However, when increasing the time interval for the change between the two 
Tsyganenko levels to 600 s (instead of 60 s), the abrupt changes vanish. This is due to increasing 
the time interval that will make the curve much smoother and the transition less abrupt. 


3.5 Dynamic Magnetic and Electric Fields 

The comparison concerning a dynamic magnetic field showed a very good resemblance between 
the two cases. When trying to add a dynamic electric field, the mathematics at first posed a 
problem. The comparison did not show enough resemblance between the two cases. The problem 
was finally solved and when comparing the results of a dynamic magnetic and electric field, the 
results of the precalculated case showed a very good resemblance with the analytical results. 
Satisfactory results were, therefore, achieved regarding a dynamic magnetic and electric field. 
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4. Conclusions 


The results from the particle motion in an interpolated magnetic field and particle motion in 
interpolated magnetic and electric fields showed a very good resemblance; therefore, we can be 
confident when it comes to executing MHD- simulations (with respect to interpolating 
simulations). A Olot of time can be saved by using interpolation — instead of time consuming 
analytical simulations — in these cases. 

The sensitivity analysis — for the particle motion in interpolated magnetic field and particle 
motions in interpolated magnetic and electric fields — in general show large improvements in 
accuracy between the analytical results and the precalculated results with an increased number of 
grid points. In some cases, more or less identical results between the analytical and the 
precalculated simulations were found. The decreased accuracy concerning energy could be 
explained with the random nature of the particle trajectories. Weighing the simulation time (the 
simulations take longer when increasing the number of grids) and the resulting accuracy, an 
optimal number of grids would probably be eight times the original number. 

The few things that differed when comparing the polynomial results with the linear results — 
dynamic magnetic and electric fields — were generally easily corrected by increasing the time 
interval between the change of the Tsyganenko levels. By doing so, the time interval will make 
the curve much smoother and the transition less abrupt. 

The comparison concerning a dynamic magnetic field showed a very good resemblance between 
the two cases, but when trying to add a dynamic electric field, the mathematics posed a problem. 
The comparison did not show enough resemblance between the two cases. The problem was how 
to calculate the expression given by 


J (E. • B)ds , (16) 

without field line tracing. Equation (16) is found in the expression for the total Zf- field, given by 

E = -V(P + J(E i »B)ds) + E u , (17) 

where P is the electrostatic potential (obtained from mapping to the ionosphere), E, is the 
inductive electric field and B is the magnetic field. 

The problem was solved in the following way: the total E'-field was precalculated during the 
transition between the two levels, calculating the E - field at a 3 second interval. The particle 
trajectory was calculated by doing temporal and spatial interpolation to find the E - field at a 
particular time and location along the trajectory. When comparing the results of a dynamic 
magnetic and electric field, a very good resemblance with the analytical results was found; thus, 
satisfactory results were achieved regarding a dynamic magnetic and electric field. 
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Figure Al. Polynomial transition, Scenario 8. 
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Figure A2. Linear transition, Scenario 8. 
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Figure A3. Polynomial transition, Scenario 8. Increasing the time interval to 600. 
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Figure A4. Linear transition, Scenario 8. Increasing the time interval to 600. 
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Figure A6. Precalculated case for Scenario 10. 
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Figure A7. Analytical case (Scenario 3). 
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Figure A8. Precalculated case with original grid numbers (Scenario 3). 
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Figure A9. Precalculated case with 8 times the number of grids (Scenario 3). 
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Figure A10. The precalculated case with 27 times the number of grids (Scenario 3). 
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Figure All. Analytical case (Scenario 16). 
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Figure A12. Precalculated case with original grid numbers (Scenario 16). 
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Figure A13. Precalculated case with 8 times the number of grids (Scenario 16). 
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Figure A14. Precalculated case with 27 times the number of grids (Scenario 16). 
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